package fem2.element;

import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IMatrixReadOnly;
import inf.jlinalg.IVector;
import inf.jlinalg.IVectorReadOnly;

import java.util.ArrayList;

import math2.LinAlgUtilities;
import fem2.Block;
import fem2.Element;
import fem2.Face;
import fem2.FieldInfo;
import fem2.IntegrationManager;
import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.NotImplementedError;
import fem2.enu.IntegrationType;
import fem2.enu.MassMatrixType;
import fem2.enu.MeshPartType;
import fem2.enu.StressType;

/**
 * 
 * Linear kinematics structural element.<br>
 * Currently, this element support linear elastic constitutive law and various
 * integration type. It wrap up all standard finite element formulation in small
 * strain regime
 * 
 * @author hbui
 * 
 */
public class StructuralElement extends Element {

	private MaterialModel mm;
	protected ArrayList<IVector> sigma = new ArrayList<IVector>();

	public StructuralElement(MeshPart mp, MaterialModel mm) {
		super(mp);
		setMaterialModel(mm);
		/*
		 * set default integration type to this element
		 */
		setIntegrationType(IntegrationType.FULL_INTEGRATION);
	}

	/**
	 * @return the face if this element associate with a face
	 */
	public Face getFace() {
		if (getMeshPart() instanceof Face)
			return (Face) getMeshPart();
		else
			throw new Error("associated mesh part was not face");
	}

	/**
	 * 
	 * @return the block if this element associate with the block
	 */
	public Block getBlock() {
		if (getMeshPart() instanceof Block)
			return (Block) getMeshPart();
		else
			throw new Error("associated mesh part was not block");
	}

	@Override
	public void setIntegrationType(IntegrationType integrationType) {
		super.setIntegrationType(integrationType);

		MeshPartType mpType = getMeshPart().getType();

		/*
		 * handling some exception for integration type
		 */
		if ((mpType == MeshPartType.TRI3 || mpType == MeshPartType.TET4)
				&& (integrationType == IntegrationType.REDUCED_INTEGRATION || integrationType == IntegrationType.SELECTIVE_REDUCED_INTEGRATION)) {
			throw new Error(
					"3-node triangle and 4-node tetrahedra don't support reduced or selective reduced integration");
		}
	}

	@Override
	public int getFieldId(int idx) {
		if (idx == 0) {
			return FieldInfo.U1;
		} else if (idx == 1) {
			return FieldInfo.U2;
		} else if (idx == 2) {
			return FieldInfo.U3;
		} else {
			throw new Error("Invalid input");
		}
	}

	/**
	 * Returns the constitutive model associated with the element.
	 * 
	 * @return the constitutive model
	 */
	public MaterialModel getMaterialModel() {
		return this.mm;
	}

	/**
	 * set the material model
	 * 
	 * @param mm
	 */
	public void setMaterialModel(MaterialModel mm) {
		this.mm = mm;
		if (mm != null) {
			mm.registerElement(this);
		}
	}

	/**
	 * Computes the standard B operator
	 * 
	 * For 2d:
	 * 
	 * <pre>
	 * 
	 * 	|dN1/dX1	0        	dN2/dX1		0	.	.	.	|
	 * B =	|0		dN1/dX2		0		dN2/dX2	.	.	.	|
	 * 	|dN1/dX2  	dN1/dX1  	dN2/dX2		dN2/dX1	.	.	.	|
	 * </pre>
	 * 
	 */
	public void computeB(MeshPart mp, IMatrix b, double... xi) {
		int dim = xi.length;
		int n = mp.getBasis().getDimension();
		double[][] gGradients = mp.gGradientsAt(xi);
		if (dim == 2) {
			for (int i = 0; i < n; i++) {
				double g1 = gGradients[i][0];
				double g2 = gGradients[i][1];

				b.set(0, 2 * i, g1);
				b.set(1, 2 * i + 1, g2);
				b.set(2, 2 * i, g2);
				b.set(2, 2 * i + 1, g1);
			}
		} else if (dim == 3) {
			for (int i = 0; i < n; i++) {
				double g1 = gGradients[i][0];
				double g2 = gGradients[i][1];
				double g3 = gGradients[i][2];

				b.set(0, 3 * i, g1);
				b.set(1, 3 * i + 1, g2);
				b.set(2, 3 * i + 2, g3);
				b.set(3, 3 * i, g2);
				b.set(3, 3 * i + 1, g1);
				b.set(4, 3 * i + 1, g3);
				b.set(4, 3 * i + 2, g2);
				b.set(5, 3 * i, g3);
				b.set(5, 3 * i + 2, g1);
			}
		}
	}

	/**
	 * Computes the standard B operator
	 * 
	 * For 2d:
	 * 
	 * <pre>
	 * 
	 *     		| dN1/dX1  0        			dN2/dX1		0       		...		|
	 * B = 	| 0      			dN1/dX2    0	    			dN2/dX2 ...		|
	 *     		| dN1/dX2  dN1/dX1  	dN2/dX2  	dN2/dX1 ...		|
	 * </pre>
	 * 
	 * */
	public static double[][] computeB(MeshPart mp, double[][] b, double... xi) {
		int dim = xi.length;
		int n = mp.getBasis().getDimension();
		double[][] gGradients = mp.gGradientsAt(xi);
		if (dim == 2) {
			for (int i = 0; i < n; i++) {
				double g1 = gGradients[i][0];
				double g2 = gGradients[i][1];

				b[0][2 * i] = g1;
				b[1][2 * i + 1] = g2;
				b[2][2 * i] = g2;
				b[2][2 * i + 1] = g1;
			}
		} else if (dim == 3) {
			for (int i = 0; i < n; i++) {
				double g1 = gGradients[i][0];
				double g2 = gGradients[i][1];
				double g3 = gGradients[i][2];

				b[0][3 * i] = g1;
				b[1][3 * i + 1] = g2;
				b[2][3 * i + 2] = g3;
				b[3][3 * i] = g2;
				b[3][3 * i + 1] = g1;
				b[4][3 * i + 1] = g3;
				b[4][3 * i + 2] = g2;
				b[5][3 * i] = g3;
				b[5][3 * i + 2] = g1;
			}
		}
		return b;
	}

	@Override
	public IMatrixReadOnly computeKte(IMatrix ke) {
		MeshPart mp = getMeshPart();
		int n = mp.getBasis().getDimension();
		int dim = mp.getDimension();
		int ntrs = dim * (dim + 1) / 2;
		Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * n);
		double th = getMaterialModel().getSectionProperty();

		if (getIntegrationType() != IntegrationType.SELECTIVE_REDUCED_INTEGRATION) {
			int ng = getNumIntegrationPoints();
			for (int k = 0; k < ng; k++) {
				double w = getIntegrationWeight(k);
				double[] xi = getIntegrationPoint(k);

				computeB(mp, B, xi);
				double[][] jac = mp.jacXAt(xi);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = w * detJ * th;

				IMatrix C = mm.getConstitutiveMatrix(this, k);
				LinAlgUtilities.btCB(B, C, coeff, ke);
			}
		} else if (getIntegrationType() == IntegrationType.SELECTIVE_REDUCED_INTEGRATION) {
			int ng = getNumIntegrationPoints();
			int ngFull = IntegrationManager.getInstance().getNumIntegrationPoints(mp.getType(),
					getIntegrationOrder());
			IMatrix tmp = new Array2DMatrix(ntrs, ntrs);
			/*
			 * compute first part using full integration
			 */
			for (int k = 0; k < ngFull; k++) {
				double w = getIntegrationWeight(k);
				double[] xi = getIntegrationPoint(k);

				computeB(mp, B, xi);
				double[][] jac = mp.jacXAt(xi);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = w * detJ * th;

				IMatrix C = mm.getConstitutiveMatrix(this, k);
				BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, IntegrationManager.getInstance()
						.getVolumeIntegralIminusDMatrix(dim), BLAM.NO_TRANSPOSE, C, 0.0, tmp);
				LinAlgUtilities.btCB(B, tmp, coeff, ke);
			}
			/*
			 * compute second part using reduced integration
			 */
			for (int k = ngFull; k < ng; k++) {
				double wReduced = getIntegrationWeight(k);
				double[] xiReduced = getIntegrationPoint(k);

				computeB(mp, B, xiReduced);
				double[][] jac = mp.jacXAt(xiReduced);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = wReduced * detJ * th;

				IMatrix C = mm.getConstitutiveMatrix(this, ng + k);
				BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, IntegrationManager.getInstance()
						.getVolumeIntegralDMatrix(dim), BLAM.NO_TRANSPOSE, C, 0.0, tmp);
				LinAlgUtilities.btCB(B, tmp, coeff, ke);
			}
		}
		return ke;
	}

	@Override
	public IVectorReadOnly computeRie(IVector rie) {
		MeshPart mp = getMeshPart();
		int n = mp.getBasis().getDimension();
		int dim = mp.getDimension();
		int ntrs = dim * (dim + 1) / 2;
		Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * n);
		double th = getMaterialModel().getSectionProperty();

		if (getIntegrationType() != IntegrationType.SELECTIVE_REDUCED_INTEGRATION) {
			int ng = getNumIntegrationPoints();
			for (int k = 0; k < ng; k++) {
				double w = getIntegrationWeight(k);
				double[] xi = getIntegrationPoint(k);

				computeB(mp, B, xi);
				double[][] jac = mp.jacXAt(xi);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = w * detJ * th;

				IVector s = getStress(k);
				BLAM.multiply(coeff, BLAM.TRANSPOSE, B, s, 1.0, rie);
			}
		} else if (getIntegrationType() == IntegrationType.SELECTIVE_REDUCED_INTEGRATION) {
			int ng = getNumIntegrationPoints();
			int ngFull = IntegrationManager.getInstance().getNumIntegrationPoints(mp.getType(),
					getIntegrationOrder());
			IVector tmp = new ArrayVector(ntrs);
			/*
			 * compute first part by full integration scheme
			 */
			for (int k = 0; k < ngFull; k++) {
				double w = getIntegrationWeight(k);
				double[] xi = getIntegrationPoint(k);

				computeB(mp, B, xi);
				double[][] jac = mp.jacXAt(xi);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = w * detJ * th;

				IVector s = getStress(k);
				BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, IntegrationManager.getInstance()
						.getVolumeIntegralIminusDMatrix(dim), s, 0.0, tmp);
				BLAM.multiply(coeff, BLAM.TRANSPOSE, B, tmp, 1.0, rie);
			}
			/*
			 * compute second part by reduced integration scheme
			 */
			for (int k = ngFull; k < ng; k++) {
				double wReduced = getIntegrationWeight(k);
				double[] xiReduced = getIntegrationPoint(k);

				computeB(mp, B, xiReduced);
				double[][] jac = mp.jacXAt(xiReduced);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = wReduced * detJ * th;

				IVector s = getStress(k);
				BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, IntegrationManager.getInstance()
						.getVolumeIntegralDMatrix(dim), s, 0.0, tmp);
				BLAM.multiply(coeff, BLAM.TRANSPOSE, B, tmp, 1.0, rie);
			}
		}
		return rie;
	}

	@Override
	public IVectorReadOnly computeRes(IVector res) {
		MeshPart mp = getMeshPart();
		int n = mp.getBasis().getDimension();
		int dim = mp.getDimension();
		int ntrs = dim * (dim + 1) / 2;
		Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * n);
		double th = getMaterialModel().getSectionProperty();

		if (getIntegrationType() != IntegrationType.SELECTIVE_REDUCED_INTEGRATION) {
			int ng = getNumIntegrationPoints();
			for (int k = 0; k < ng; k++) {
				double w = getIntegrationWeight(k);
				double[] xi = getIntegrationPoint(k);

				computeB(mp, B, xi);
				double[][] jac = mp.jacXAt(xi);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = w * detJ * th;

				IVector s = getStress(k, StressType.CAUCHY_STRESS);
				BLAM.multiply(coeff, BLAM.TRANSPOSE, B, s, 1.0, res);
			}
		} else if (getIntegrationType() == IntegrationType.SELECTIVE_REDUCED_INTEGRATION) {
			int ng = getNumIntegrationPoints();
			int ngFull = IntegrationManager.getInstance().getNumIntegrationPoints(mp.getType(),
					getIntegrationOrder());
			IVector tmp = new ArrayVector(ntrs);
			/*
			 * compute first part by full integration scheme
			 */
			for (int k = 0; k < ngFull; k++) {
				double w = getIntegrationWeight(k);
				double[] xi = getIntegrationPoint(k);

				computeB(mp, B, xi);
				double[][] jac = mp.jacXAt(xi);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = w * detJ * th;

				IVector s = getStress(k, StressType.CAUCHY_STRESS);
				BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, IntegrationManager.getInstance()
						.getVolumeIntegralIminusDMatrix(dim), s, 0.0, tmp);
				BLAM.multiply(coeff, BLAM.TRANSPOSE, B, tmp, 1.0, res);
			}
			/*
			 * compute second part by reduced integration scheme
			 */
			for (int k = ngFull; k < ng; k++) {
				double wReduced = getIntegrationWeight(k);
				double[] xiReduced = getIntegrationPoint(k);

				computeB(mp, B, xiReduced);
				double[][] jac = mp.jacXAt(xiReduced);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = wReduced * detJ * th;

				IVector s = getStress(k, StressType.CAUCHY_STRESS);
				BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, IntegrationManager.getInstance()
						.getVolumeIntegralDMatrix(dim), s, 0.0, tmp);
				BLAM.multiply(coeff, BLAM.TRANSPOSE, B, tmp, 1.0, res);
			}
		}
		return res;
	}

	@Override
	public IMatrixReadOnly computeMe(IMatrix me, MassMatrixType massType) {
		IntegrationManager integrator = IntegrationManager.getInstance();
		MeshPart mp = getMeshPart();
		int n = mp.getBasis().getDimension();
		int dim = mp.getDimension();
		double th = getMaterialModel().getSectionProperty();
		double rho = getMaterialModel().getDensity();
		double s = 1.0;

		if (getIntegrationType() == IntegrationType.SELECTIVE_REDUCED_INTEGRATION) {
			s = 0.5; // to account for more Gauss points in the bulk
		}

		if (massType == MassMatrixType.VARIATIONAL_MASS_LUMPING) {
			/*
			 * this variational mass lumping used a stiffness-consistent
			 * approach
			 */

			/*
			 * use integration scheme one higher order than integration scheme
			 * for stiffness
			 */
			int p1 = getIntegrationOrder() + 1;
			int ng = integrator.getNumIntegrationPoints(mp.getType(), p1);
			for (int k = 0; k < ng; k++) {
				double w = integrator.getIntegrationWeight(mp.getType(), p1, k);
				double[] xi = integrator.getIntegrationPoint(mp.getType(), p1, k);

				double[][] jac = mp.jacXAt(xi);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = w * detJ * th * rho * s;

				double[] N = mp.valueAt(xi);

				for (int i = 0; i < n; i++) {
					for (int j = 0; j < n; j++) {
						for (int l = 0; l < dim; l++) {
							me.add(i * dim + l, j * dim + l, coeff * N[i] * N[j]);
						}
					}
				}
			}
		} else if (massType == MassMatrixType.DIRECT_MASS_LUMPING) {
			double m = rho * mp.computeVolume() * th / n;
			for (int i = 0; i < n; i++) {
				for (int j = 0; j < dim; j++) {
					int k = i * dim + j;
					me.set(k, k, m);
				}
			}
		} else {
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		}

		return me;
	}

	/**
	 * compute geometrical stiffness matrix
	 * 
	 * @param ke
	 * @return
	 */
	public IMatrixReadOnly computeKge(IMatrix ke) {
		MeshPart mp = getMeshPart();
		int n = mp.getBasis().getDimension();
		int dim = mp.getDimension();
		int ng = getNumIntegrationPoints();
		double th = getMaterialModel().getSectionProperty();
		double s = 1.0;

		if (getIntegrationType() == IntegrationType.SELECTIVE_REDUCED_INTEGRATION) {
			s = 0.5; // to account for more Gauss points in the bulk
		}

		for (int i = 0; i < ng; i++) {
			double[] xi = getIntegrationPoint(i);
			double w = getIntegrationWeight(i);

			double[][] jac = mp.jacXAt(xi);
			double detJ = LinAlgUtilities.determinantNxN(jac);
			double[][] gv = mp.gGradientsAt(xi);
			IVector sigma = getStress(i);

			if (dim == 2) {
				for (int j = 0; j < n; j++) {
					for (int k = 0; k < n; k++) {

						double ni1 = gv[j][0];
						double ni2 = gv[j][1];
						double nj1 = gv[k][0];
						double nj2 = gv[k][1];

						double f = sigma.get(0) * ni1 * nj1 + sigma.get(1) * ni2 * nj2
								+ sigma.get(2) * (ni2 * nj1 + ni1 * nj2);
						double coeff = w * detJ * f * th * s;
						ke.add(2 * j, 2 * k, coeff);
						ke.add(2 * j + 1, 2 * k + 1, coeff);
					}
				}
			} else if (dim == 3) {
				for (int j = 0; j < n; j++) {
					for (int k = 0; k < n; k++) {

						double ni1 = gv[j][0];
						double ni2 = gv[j][1];
						double ni3 = gv[j][2];
						double nj1 = gv[k][0];
						double nj2 = gv[k][1];
						double nj3 = gv[k][2];

						double f = sigma.get(0) * ni1 * nj1 + sigma.get(1) * ni2 * nj2
								+ sigma.get(2) * ni3 * nj3 + sigma.get(3) * (ni2 * nj1 + ni1 * nj2)
								+ sigma.get(4) * (ni2 * nj3 + ni3 * nj2) + sigma.get(5)
								* (ni3 * nj1 + ni1 * nj3);
						double coeff = w * detJ * f * th * s;
						ke.add(3 * j, 3 * k, coeff);
						ke.add(3 * j + 1, 3 * k + 1, coeff);
						ke.add(3 * j + 2, 3 * k + 2, coeff);
					}
				}
			}
		}
		return ke;
	}

	/**
	 * 
	 * @param xi
	 * @return regular strain at natural coordinates xi by B*ue
	 */
	public IVector computeRegularStrain(double... xi) {
		MeshPart mp = getMeshPart();
		IVector u = new ArrayVector(mp.getUHatLocal());
		int dim = mp.getDimension();
		int n = mp.getBasis().getDimension();
		int ntrs = dim * (dim + 1) / 2;
		Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * n);
		computeB(mp, B, xi);
		IVector e = new ArrayVector(ntrs);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, B, u, 0.0, e);
		return e;
	}

	/**
	 * 
	 * Computes the strain at natural coordinate xi.
	 * 
	 * @param xi
	 * @return
	 */
	public IVector computeStrain(double... xi) {
		return computeRegularStrain(xi);
	}

	// /**
	// * Computes the strain at Gauss point k.
	// *
	// * @param k
	// * @return
	// */
	// public IVector computeStrain(int k) {
	// return computeStrain(getGaussPoint(k));
	// }

	/**
	 * get the (algorithmic) stress at Gauss point k
	 * 
	 * @param k
	 * @return the sigma
	 */
	public IVector getStress(int k) {
		return sigma.get(k);
	}

	/**
	 * get the stress at Gauss point k specified by type (algorithmic, recover,
	 * ...)
	 * 
	 * @param k
	 * @return the sigma
	 */
	public IVector getStress(int k, StressType stressType) {
		switch (stressType) {
		case ALGORITHMIC_STRESS:
			return getStress(k);
		default:
			return getMaterialModel().getStress(this, k, stressType);
		}
	}

	@Override
	public int getNumIntegrationPoints() {
		switch (getIntegrationType()) {
		case FULL_INTEGRATION:
			return IntegrationManager.getInstance().getNumIntegrationPoints(
					getMeshPart().getType(), getIntegrationOrder());
		case REDUCED_INTEGRATION:
			return IntegrationManager.getInstance().getNumIntegrationPoints(
					getMeshPart().getType(), getIntegrationOrder() - 1);
		case SELECTIVE_REDUCED_INTEGRATION:
			return IntegrationManager.getInstance().getNumIntegrationPoints(
					getMeshPart().getType(), getIntegrationOrder())
					+ IntegrationManager.getInstance().getNumIntegrationPoints(
							getMeshPart().getType(), getIntegrationOrder() - 1);
		default:
			return IntegrationManager.getInstance().getNumIntegrationPoints(this);
		}
	}

	@Override
	public double[] getIntegrationPoint(int k) {
		switch (getIntegrationType()) {
		case FULL_INTEGRATION:
			return IntegrationManager.getInstance().getIntegrationPoint(getMeshPart().getType(),
					getIntegrationOrder(), k);
		case REDUCED_INTEGRATION:
			return IntegrationManager.getInstance().getIntegrationPoint(getMeshPart().getType(),
					getIntegrationOrder() - 1, k);
		case SELECTIVE_REDUCED_INTEGRATION:
			int ng = IntegrationManager.getInstance().getNumIntegrationPoints(
					getMeshPart().getType(), getIntegrationOrder());
			if (k < ng) {
				/*
				 * return Gauss point at the full part
				 */
				return IntegrationManager.getInstance().getIntegrationPoint(
						getMeshPart().getType(), getIntegrationOrder(), k);
			} else {
				/*
				 * return Gauss point at the reduced part
				 */
				return IntegrationManager.getInstance().getIntegrationPoint(
						getMeshPart().getType(), getIntegrationOrder() - 1, k - ng);
			}
		default:
			return IntegrationManager.getInstance().getIntegrationPoint(this, k);
		}
	}

	@Override
	public double getIntegrationWeight(int k) {
		switch (getIntegrationType()) {
		case FULL_INTEGRATION:
			return IntegrationManager.getInstance().getIntegrationWeight(getMeshPart().getType(),
					getIntegrationOrder(), k);
		case REDUCED_INTEGRATION:
			return IntegrationManager.getInstance().getIntegrationWeight(getMeshPart().getType(),
					getIntegrationOrder() - 1, k);
		case SELECTIVE_REDUCED_INTEGRATION:
			int ng = IntegrationManager.getInstance().getNumIntegrationPoints(
					getMeshPart().getType(), getIntegrationOrder());
			if (k < ng) {
				/*
				 * return Gauss weight at the full part
				 */
				return IntegrationManager.getInstance().getIntegrationWeight(
						getMeshPart().getType(), getIntegrationOrder(), k);
			} else {
				/*
				 * return Gauss weight at the reduced part
				 */
				return IntegrationManager.getInstance().getIntegrationWeight(
						getMeshPart().getType(), getIntegrationOrder() - 1, k - ng);
			}
		default:
			return IntegrationManager.getInstance().getIntegrationWeight(this, k);
		}
	}

	@Override
	public void Initialize() {
		/*
		 * initialize the sigma array for this element
		 */
		int ng = getNumIntegrationPoints();
		int dim = getMeshPart().getDimension();
		sigma.clear();
		for (int i = 0; i < ng; i++) {
			sigma.add(new ArrayVector(dim * (dim + 1) / 2));
		}
		/*
		 * initialize material
		 */
		mm.Initialize(this);
	}

	@Override
	public void InitializeSolutionStep() {
		mm.InitializeSolutionStep(this);
	}

	@Override
	public void InitalizeNonlinearIteration() {
		mm.InitalizeNonlinearIteration(this);
	}

	@Override
	public void FinalizeNonlinearIteration() {
		mm.FinalizeNonlinearIteration(this);
	}

	@Override
	public void FinalizeSolutionStep() {
		mm.FinalizeSolutionStep(this);
	}

	@Override
	public void CalculateOutputData() {
	}

}
